High-Accuracy Neural Network Interatomic Potential for Silicon Nitride

In the field of machine learning (ML) and data science, it is meaningful to use the advantages of ML to create reliable interatomic potentials. Deep potential molecular dynamics (DEEPMD) are one of the most useful methods to create interatomic potentials. Among ceramic materials, amorphous silicon nitride (SiNx) features good electrical insulation, abrasion resistance, and mechanical strength, which is widely applied in industries. In our work, a neural network potential (NNP) for SiNx was created based on DEEPMD, and the NNP is confirmed to be applicable to the SiNx model. The tensile tests were simulated to compare the mechanical properties of SiNx with different compositions based on the molecular dynamic method coupled with NNP. Among these SiNx, Si3N4 has the largest elastic modulus (E) and yield stress (σs), showing the desired mechanical strength owing to the largest coordination numbers (CN) and radial distribution function (RDF). The RDFs and CNs decrease with the increase of x; meanwhile, E and σs of SiNx decrease when the proportion of Si increases. It can be concluded that the ratio of nitrogen to silicon can reflect the RDFs and CNs in micro level and macro mechanical properties of SiNx to a large extent.


Introduction
Molecular dynamics (MD) mainly rely on Newtonian mechanics to simulate the motion of molecular systems. Compared with the expensive cost of the experiment, the MD method is a cost-effective and efficient tool for exploring the properties of various complicated new materials. To make sure the simulation results match well with the experiments, it is important to use an accurate description of interatomic interactions. Quantum mechanics (QM) simulations, such as the ab initio molecular dynamics (AIMD) method based on the density functional theory (DFT), are the most reliable way to describe the atomic interactions for different systems [1][2][3][4][5]. Although AIMD exhibits desirable computational accuracy, the time cost of AIMD is very high, which limits the application of AIMD. To balance the calculation performance and speed of MD simulations, empirical interatomic potentials have been applied to MD simulations. Empirical interatomic potentials, including Lennard-Jones (LJ) [6], embedded atom method (EAM) [7,8], the Stillinger-Weber (SW) [9,10], Tersoff [11] and charge-optimized many-body (COMB) [12] potentials obviously improved the calculation performance and speed of MD. However, the application of these empirical potentials is hindered by their poor transferability. In terms of systems described by two-body interactions, the LJ potential has favorable accuracy. The SW and Tersoff potentials are able to combine two-body and three-body interactions to stabilize tetrahedral solids, but the description accuracy of bond breaking and metallic phases of silicon and carbon is not sufficient. Owing to the rapid development of machine learning (ML) and data science, it is meaningful to take the advantages of ML to create reliable ML interatomic potentials to replace the conventional empirical potentials [13][14][15]. In 2007,

Materials and Methods
DEEPMD is one of the most popular methods to create interatomic potentials recently owing to its calculation speed and accuracy [17,33]. The details of the NNP training process are listed as follows.

AIMD Calculations
In our work, Si 3 N 4 was chosen as an example to train the NNP because of the high stability of the 3:4 composition. The training database of Si 3 N 4 at different temperatures was obtained by first-principles calculations using the Vienna ab initio simulation package (VASP) [34]. The exchange-correlation interaction was described by Perdew-Burke-Ernzerhof (PBE) functional [35]. The interaction between electrons and ions was described by the projector-augmented wave (PAW) approach. The cut-off energy is 520 eV. We use the k-point mesh grid with a spacing of 0.4 Å within the Gamma-centered k-sampling to sample the Brillouin zone. The initial Si 3 N 4 configuration with 112 atoms is a cube with randomly distributed atoms; the model was built by LAMMPS "create_atoms" command and modified by Material Studio. The periodic boundary conditions were applied in x, y, and z directions. The cube size in the x, y, and z directions is 11.935 Å, 11.935 Å, and 11.55 Å, respectively. The AIMD calculations were fulfilled at a constant temperature of 2000 K with an NVT ensemble. The Nose-Hoover thermostat was used to control the temperature of the AIMD simulation. The timestep is 1 fs running 10,000 steps. The energy and force errors less than 10 −5 eV/atom and 0.01 eV/Å, respectively, are convergence criteria for geometry optimization.

Deep Potential Training Process
Based on the AIMD calculations, 10,000 data points were transferred from the output file "OUTCAR"; the data points were divided into five sets. The four sets were used for training databases, with the remaining one selected as a testing database. The smooth edition of the DEEPMD and DeepPot-SE model, implemented in the DEEPMD-kit package [33], was used to train the interatomic potential. The cutoff radius of the model is 6.0 Å for neighbor searching with the smoothing function starting from 5.8 Å. The hidden layers were divided into three layers, and the number of neurons in each layer is 25, 50, and 100, respectively. The learning rate starts from 0.001 using a decay rate of 0.95 every 5000 steps. The decay step and stop learning rate are 5000 and 3.51 × 10 −8 , respectively. Based on the datasets of AIMD calculations at 2000 K, the NNP potential at 2000 K was obtained. In order to obtain a high-quality potential, a training database including different temperatures is crucial; as a result, the same method was used for the 3000 K potential training as well. Then, the 2000 K and 3000 K NNP potentials were combined to achieve a new NNP. Finally, the NNP was frozen, and the frozen model can be used in model testing and MD simulations.

The Accuracy of NNP
To verify the accuracy of NNP at 2000 K and 3000 K, the energies and forces from NNP after DEEPMD training were compared with those from AIMD shown in Figures 1 and 2. The results show that the validation data generally distribute around the y = x line, showing that the obtained potential can predict the energies and forces precisely. The root-mean-square errors (RMSEs) and R-Square (R 2 ) of energies were calculated to evaluate the performance of NNP. The smaller RMSEs mean better prediction. On the contrary, the higher R 2 , between 0 and 1, indicates better results. It can be seen that the NNP with larger R 2 and smaller RMSE at y and z directions, respectively, works better than that of the x direction, both in Figures 1 and 2. Meanwhile, the R 2 and RMSE are acceptable for NNP at x direction as well.
To further validate the reliability of NNP, the RDFs, CNs, and BADs of Si 3 N 4 calculated by NNP were compared with those of AIMD. To research these properties of Si 3 N 4 , a simulation at 2000 K was conducted by a large-scale atomic/molecular massively parallel simulator (LAMMPS) [36] with the MD [37]. The detailed simulation parameters are shown in the Supplementary Material. The RDFs, also known as pair correlation function, usually refer to the distribution probability of other particles in the ∆r thickness shell at the distance r of a specified particle. RDFs are widely used to study the degree of order of materials and describe the correlation of atoms, which can be calculated by where V is the volume of the simulation cell. N α and N β are the number of α-type ions and β-type ions, respectively. N αβ (r, ∆r) is the average number of α-type ions around β-type ions in a spherical space. The results of the NNP + MD and AIMD simulations in Figure 3a confirm the accuracy of the NNP. The MD in the following content is specified as MD coupled with NNP. It can be seen that the first peak values of Si 3 N 4 RDFs calculated by MD are located around r = 1.7 Å. The results show that the distribution function g(r) obtained by MD is consistent with the AIMD results calculated by VASP. The same conclusion can be reached at 3000 K in Figure 4a. Therefore, it is concluded that the NNP is capable of predicting structure information of SiN x with AIMD accuracy within this range from 2000 K to 3000 K.       coupled with NNP. It can be seen that the first peak values of Si3N4 RDFs calculated by MD are located around r = 1.7 Å. The results show that the distribution function g(r) obtained by MD is consistent with the AIMD results calculated by VASP. The same conclusion can be reached at 3000 K in Figure 4a. Therefore, it is concluded that the NNP is capable of predicting structure information of SiNx with AIMD accuracy within this range from 2000 K to 3000 K. ions in a spherical space. The results of the NNP + MD and AIMD simulations in Figure  3a confirm the accuracy of the NNP. The MD in the following content is specified as MD coupled with NNP. It can be seen that the first peak values of Si3N4 RDFs calculated by MD are located around r = 1.7 Å. The results show that the distribution function g(r) obtained by MD is consistent with the AIMD results calculated by VASP. The same conclusion can be reached at 3000 K in Figure 4a. Therefore, it is concluded that the NNP is capable of predicting structure information of SiNx with AIMD accuracy within this range from 2000 K to 3000 K. CNs are the number of coordination atoms around the central atom of a compound. The RDFs depend on the multilayered coordination radius and coordination particle numbers. The CNs could be calculated by integration of the RDFs where N is the CN. ρ is particle density in the simulation cell. R min is the first minimum in Figures 3a and 4a. As is shown in Figure 3b, we now concentrate on the CNs of Si 3 N 4 . The results calculated by MD are compared with the AIMD simulation results. The CNs calculated by MD fit the AIMD results accurately, as well as in Figure 4b. The CNs are proportional to the pair distance and CNs are zero when the pair distance is less than 1.5 Å, so the minimal pair distance of Si 3 N 4 is 1.5 Å. The BADs were calculated to analyze the local geometries of the first coordination shell, which can be calculated by where θ is the bond angle. r ij , r ik, and r jk are the bond length between atoms. The BADs of Si 3 N 4 obtained by MD and AIMD at 2000 K and 3000 K are shown in Figures 5 and 6 localized between 75 • and 120 • , and the Si-N-Si and N-Si-N peak values of MD BADs are 90 • , which is consistent with AIMD BADs. As for 3000 K, though the curves are rougher than that of 2000 K, the Si-N-Si and N-Si-N peak values of MD BADs are consistent well with AIMD BADs. Therefore, NNP is reliable enough to predict CNs and BADs information of SiN x . Based on the RDFs, CNs, and BADs analysis above, it can be concluded that we can use the NNP to calculate the structure information of SiN x .
local geometries of the first coordination shell, which can be calculated by where θ is the bond angle. rij, rik, and rjk are the bond length between atoms. The BADs of Si3N4 obtained by MD and AIMD at 2000 K and 3000 K are shown in Figure 5 and Figure  6, respectively. At 2000 K, there are few BAD smaller than 45°. The majority of BADs are localized between 75° and 120°, and the Si-N-Si and N-Si-N peak values of MD BADs are 90°, which is consistent with AIMD BADs. As for 3000 K, though the curves are rougher than that of 2000 K, the Si-N-Si and N-Si-N peak values of MD BADs are consistent well with AIMD BADs. Therefore, NNP is reliable enough to predict CNs and BADs information of SiNx. Based on the RDFs, CNs, and BADs analysis above, it can be concluded that we can use the NNP to calculate the structure information of SiNx.

The NPP Applied for Structure Information
Similar to the empirical potentials, the NNP can also be used for different compositions with the same elements. To study the effect of compositions on SiNx structures, we used the NNP to simulate the heating process of SiNx and analyzed the simulation results, including RDFs, CNs, and BADs. The models of SiNx were built by LAMMPS, as shown in Figure 7. The parameters of SiNx models are shown in Table 1. The boundary condition, timestep, ensemble, and other parameters of MD simulations are the same as part 3 for Si3N4.

The NPP Applied for Structure Information
Similar to the empirical potentials, the NNP can also be used for different compositions with the same elements. To study the effect of compositions on SiN x structures, we used the NNP to simulate the heating process of SiN x and analyzed the simulation results, including RDFs, CNs, and BADs. The models of SiN x were built by LAMMPS, as shown in Figure 7. The parameters of SiN x models are shown in Table 1. The boundary condition, timestep, ensemble, and other parameters of MD simulations are the same as part 3 for Si 3 N 4 . tions with the same elements. To study the effect of compositions on SiNx structures, we used the NNP to simulate the heating process of SiNx and analyzed the simulation results, including RDFs, CNs, and BADs. The models of SiNx were built by LAMMPS, as shown in Figure 7. The parameters of SiNx models are shown in Table 1. The boundary condition, timestep, ensemble, and other parameters of MD simulations are the same as part 3 for Si3N4.  The RDFs of SiNx at 2000 K and 3000 K are shown in Figure 8. In terms of a given temperature, the g(r) peak values decrease while the proportion of Si composition increase, indicating that the pair distance decrease in SiNx with a higher Si proportion. In The RDFs of SiN x at 2000 K and 3000 K are shown in Figure 8. In terms of a given temperature, the g(r) peak values decrease while the proportion of Si composition increase, indicating that the pair distance decrease in SiN x with a higher Si proportion. In the case of 2000 K, the peak and valley values fluctuate more obviously. As the temperature increases to 3000 K, there are few differences between the peak values of different compositions, indicating that all the SiN x are in the same phase and there is less difference in microstructure as the temperature increases. The heating process for SiN x was calculated as well with the RDFs at different temperatures shown in Figure S1. the case of 2000 K, the peak and valley values fluctuate more obviously. As the temperature increases to 3000 K, there are few differences between the peak values of different compositions, indicating that all the SiNx are in the same phase and there is less difference in microstructure as the temperature increases. The heating process for SiNx was calculated as well with the RDFs at different temperatures shown in Figure S1. During the heating process, CNs and BADs at 3000 K are shown in Figure 9. The CNs of SiNx in Figure 9a share the same profile, which is all proportional to the pair distance. Among these CNs, SiNx, with a 3:4 composition, has the largest CN, which corresponds to the largest g(r) in the Si-poor model. Therefore, Si3N4 is expected to deliver the most stable configuration at the micro level. With the increase of x in SiNx, the CNs decrease slightly, indicating that the CNs are weakly affected by composition. The BADs of SiNx at 3000 K are shown in Figure 9b,c. Similar to the CNs, the composition has little influence on the BADs, especially for the N-Si-N BADs in Figure 9c. In summary, the RDFs, CNs, and BADs of SiNx are different, so the structures of SiNx are different as well. In the next section, we used tensile simulations to explain the impact of different compositions on SiNx mechanical properties from a macro perspective. During the heating process, CNs and BADs at 3000 K are shown in Figure 9. The CNs of SiN x in Figure 9a share the same profile, which is all proportional to the pair distance. Among these CNs, SiN x, with a 3:4 composition, has the largest CN, which corresponds to the largest g(r) in the Si-poor model. Therefore, Si 3 N 4 is expected to deliver the most stable configuration at the micro level. With the increase of x in SiN x , the CNs decrease slightly, indicating that the CNs are weakly affected by composition. The BADs of SiN x at 3000 K are shown in Figure 9b,c. Similar to the CNs, the composition has little influence on the BADs, especially for the N-Si-N BADs in Figure 9c. In summary, the RDFs, CNs, and BADs of SiN x are different, so the structures of SiN x are different as well. In the next section, we used tensile simulations to explain the impact of different compositions on SiN x mechanical properties from a macro perspective.

The NNP Applied for Tensile Tests
To evaluate the accuracy of tensile by MD, the simulated strain-stress curve of Si3N4 compared with the experimental result is shown in Figure 10. The elastic modulus E of the simulation and experiment is 284.6 GPa and 253.3 GPa, respectively, which are consistent well with the reported E in the previous reference [38,39]. As is known, the E depends on the size of the model [39,40]; the simulated E 10.9% higher than the experimental

The NNP Applied for Tensile Tests
To evaluate the accuracy of tensile by MD, the simulated strain-stress curve of Si 3 N 4 compared with the experimental result is shown in Figure 10. The elastic modulus E of the simulation and experiment is 284.6 GPa and 253.3 GPa, respectively, which are consistent well with the reported E in the previous reference [38,39]. As is known, the E depends on the size of the model [39,40]; the simulated E 10.9% higher than the experimental case is Nanomaterials 2023, 13, 1352 9 of 12 acceptable for the tensile test. It is impossible to fabricate the perfect single crystal Si 3 N 4 in the experiment. The Si 3 N 4 mechanical strength is determined by the defects and grain boundaries, so the experimental mechanical strength of Si 3 N 4 is typically smaller than those of calculated results. The tensile tests of Si 3 N 4 at different temperatures were calculated as well shown in Figure S2. The results show that the elastic modulus and peak values vary inversely with the temperature increasing, so when the temperature is 300 K, the mechanical properties of Si 3 N 4 perform better than those at 1000 K and 3000 K. Since the accuracy of MD has been confirmed, the tensile tests of SiN x with compositions from 3:4 to 1:1 were calculated at 300 K. The tensile models of SiN x were repeated to supercells from the stable structures in part 4 with parameters of tensile models shown in Figure S3. The cross-sections of SiN x in the tensile tests are shown in Figure 11. It can be seen that the cross-section of Si 3 N 4 is flat, and the strain is small, so the Si 3 N 4 cracked immediately and exhibited brittleness properties. In terms of other SiN x , the lengths in the y direction are larger than that of Si 3 N 4 when the fractures occurred, indicating that it takes a long time for the SiN x to crack. With the increase of x, the SiN x starts to exhibit flexibility, especially the SiN with the largest deformation when the fracture occurred.  Figure S2. The results show that the elastic modulus and peak values vary inversely with the temperature increasing, so when the temperature is 300 K, the mechanical properties of Si3N4 perform better than those at 1000 K and 3000 K. Since the accuracy of MD has been confirmed, the tensile tests of SiNx with compositions from 3:4 to 1:1 were calculated at 300 K. The tensile models of SiNx were repeated to supercells from the stable structures in part 4 with parameters of tensile models shown in Figure S3. The cross-sections of SiNx in the tensile tests are shown in Figure 11. It can be seen that the cross-section of Si3N4 is flat, and the strain is small, so the Si3N4 cracked immediately and exhibited brittleness properties. In terms of other SiNx, the lengths in the y direction are larger than that of Si3N4 when the fractures occurred, indicating that it takes a long time for the SiNx to crack. With the increase of x, the SiNx starts to exhibit flexibility, especially the SiN with the largest deformation when the fracture occurred.  The strain-stress curves and mechanical properties of SiNx are shown in Figure 12. All the curves have a linear increase at the initial stage, which corresponds to an elastic deformation. The slopes of SiNx at the initial linear stage are different, indicating that the E varies with different compositions. The E is an important parameter of materials at the macro level, which represents the ability of an object to resist elastic deformation and case is acceptable for the tensile test. It is impossible to fabricate the perfect single crystal Si3N4 in the experiment. The Si3N4 mechanical strength is determined by the defects and grain boundaries, so the experimental mechanical strength of Si3N4 is typically smaller than those of calculated results. The tensile tests of Si3N4 at different temperatures were calculated as well shown in Figure S2. The results show that the elastic modulus and peak values vary inversely with the temperature increasing, so when the temperature is 300 K, the mechanical properties of Si3N4 perform better than those at 1000 K and 3000 K. Since the accuracy of MD has been confirmed, the tensile tests of SiNx with compositions from 3:4 to 1:1 were calculated at 300 K. The tensile models of SiNx were repeated to supercells from the stable structures in part 4 with parameters of tensile models shown in Figure S3. The cross-sections of SiNx in the tensile tests are shown in Figure 11. It can be seen that the cross-section of Si3N4 is flat, and the strain is small, so the Si3N4 cracked immediately and exhibited brittleness properties. In terms of other SiNx, the lengths in the y direction are larger than that of Si3N4 when the fractures occurred, indicating that it takes a long time for the SiNx to crack. With the increase of x, the SiNx starts to exhibit flexibility, especially the SiN with the largest deformation when the fracture occurred.  The strain-stress curves and mechanical properties of SiNx are shown in Figure 12. All the curves have a linear increase at the initial stage, which corresponds to an elastic deformation. The slopes of SiNx at the initial linear stage are different, indicating that the E varies with different compositions. The E is an important parameter of materials at the macro level, which represents the ability of an object to resist elastic deformation and The strain-stress curves and mechanical properties of SiN x are shown in Figure 12. All the curves have a linear increase at the initial stage, which corresponds to an elastic deformation. The slopes of SiN x at the initial linear stage are different, indicating that the E varies with different compositions. The E is an important parameter of materials at the macro level, which represents the ability of an object to resist elastic deformation and reflects the bond strength between atoms, ions, or molecules at the micro level. The Si 3 N 4 curve has a steep slope and large yield stress σ s in the linear stage, illustrating that Si 3 N 4 is a typical brittle material. On the contrary, another SiN x demonstrates ductile property, especially the SiN with an obvious yield stage. When the strain ranges from 10% to 22%, the SiN curve becomes flat in the yield stage. The flexibility of SiN significantly improved compared with its Si 3 N 4 counterpart, which is consistent with the analysis in Figure 8. The E of the maximal and minimal slope is 349.78 GPa and 138.39 GPa for x = 4/3 and x = 1/1, respectively. The E is 199.68 GPa, 188.95 GPa, 240.82 GPa, and 160.48 GPa, respectively, when x = 5/4, 6/5, 7/6, and 8/7. Besides the E, the σ s of SiN x are different as well, showing that fracture occurs at different stages during the tensile simulations. The maximal σ s is 27.45 GPa of Si 3 N 4 . Although the σ s of other SiN x decrease with higher x, the σ s of SiN x fail to vary inversely with the x due to the amorphous structures of SiN x (x = 5/4, 6/5, 7/6, and 8/7). Among these curves, the 3:4 composition has the maximal E and σ s , which can correspond to the largest RDFs and CNs in Figures 8 and 9 reflects the bond strength between atoms, ions, or molecules at the micro level. The Si3N4 curve has a steep slope and large yield stress σs in the linear stage, illustrating that Si3N4 is a typical brittle material. On the contrary, another SiNx demonstrates ductile property, especially the SiN with an obvious yield stage. When the strain ranges from 10% to 22%, the SiN curve becomes flat in the yield stage. The flexibility of SiN significantly improved compared with its Si3N4 counterpart, which is consistent with the analysis in Figure 8.  Figure 8 and Figure 9, respectively. The RDFs and CNs decrease with the decreasing of x, leading to the flexibility improvement of SiNx, since for Si3N4 x = 4/3 while for SiN x = 1. Meanwhile, the E of SiNx excluding Si3N4 decrease compared with that of Si3N4, which is generally inversely proportional to the RDFs.
The results show that the RDFs and CNs at the micro level can reflect the macro mechanical properties of SiNx to a large extent.

Conclusions
In this work, the interatom potential for SiNx was created by the DEEPMD kit with a neural network. Based on the comparison of energies and forces from NNP and AIMD for tranSi3N4, we find that NNP can easily achieve DFT accuracy. The RDFs, CNs, and BADs simulations between the NNP and AIMD confirmed that the NNP is reliable enough to calculate the structure information Si3N4 as well. Then, we conducted comprehensive calculations to predict the RDFs, CNs, and BADs of SiNx. Therefore, we used tensile simulations to explain the impact of different compositions on SiNx properties from a macro perspective. Si3N4 is a typical brittle material with the largest E and σs. The flexibility of SiNx, excluding Si3N4, improved, leading to the decrease of E and σs. The E and σs fail to be inversely proportional with x due to the amorphous structures of SiNx. The RDFs and CNs at the micro level can reflect the macro mechanical properties of SiNx to a large extent. Among these compositions, x = 4/3 features high mechanical strength, owing to the largest CN and RDF. The E of SiNx, excluding Si3N4, decreases compared with that of Si3N4, leading to the flexibility improvement of SiNx, which is generally inversely proportional to the RDFs.

Conclusions
In this work, the interatom potential for SiN x was created by the DEEPMD kit with a neural network. Based on the comparison of energies and forces from NNP and AIMD for tranSi 3 N 4 , we find that NNP can easily achieve DFT accuracy. The RDFs, CNs, and BADs simulations between the NNP and AIMD confirmed that the NNP is reliable enough to calculate the structure information Si 3 N 4 as well. Then, we conducted comprehensive calculations to predict the RDFs, CNs, and BADs of SiN x . Therefore, we used tensile simulations to explain the impact of different compositions on SiN x properties from a macro perspective. Si 3 N 4 is a typical brittle material with the largest E and σ s . The flexibility of SiN x, excluding Si 3 N 4, improved, leading to the decrease of E and σ s . The E and σ s fail to be inversely proportional with x due to the amorphous structures of SiN x . The RDFs and CNs at the micro level can reflect the macro mechanical properties of SiN x to a large extent. Among these compositions, x = 4/3 features high mechanical strength, owing to the largest CN and RDF. The E of SiN x, excluding Si 3 N 4, decreases compared with that of Si 3 N 4 , leading to the flexibility improvement of SiN x , which is generally inversely proportional to the RDFs.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nano13081352/s1, Figure S1: The RDFs of SiN x with different compositions in heating process; Figure S2: Amorphous Si 3 N 4 tensile simulations; Figure S3: The whole process of fracture during tensile simulations; Table S1: The specific simulation parameters in SiN x tensile simulations.